Section 1.2, Probelm 2

To make numbers easier to deal with, I shifted the start time to zero and population is in unit of thousands.

year=Table[i, {i, 1790,2010,10}]
time= year -1790
pop={3929, 5308, 7240, 9638,12866, 17069, 23192, 31443, 38558, 50156, 62948,75995, 91972, 105711, 122755, 131669, 150697, 179323, 203212, 226505, 248710, 281416, 308746}

Here I combine the data into ordered pairs.

data=Transpose[{time, pop}]

TableForm[data]

Below is the plot of the data.

dataplot=ListPlot[data]

pop[[n]] gives the nth data point in pop. Below I have formed the first differences.

rather than counting the number of data points, I ask Mathematica to give me the length of the time vector.

Length[time]

popchange=Table[pop[[n+1]]-pop[[n]],{n,1,22}]

I will add a blank spot to the end of this list so it will be the same length as the other data and I use the transpose command to display all in table format.

popchange=Append[popchange,{}]

table=Transpose[{year, time,pop,popchange}]//TableForm

It looks like p[n+1]-p[n] is mostly increasing as the population size is increasing. One way of modeling this is to say p[n+1]-p[n] is proportional to the p[n]; p[n+1]-p[n] = k p [n]. If this is the case, the ratio (p[n+1]-p[n]) / p[n] should be the constant k. Plot (p[n+1]-p[n]) versus p[n].

conjecturedata=Table[{pop[[n]],pop[[n+1]]-pop[[n]]},{n,1,22}]

TableForm[conjecturedata]

ListPlot[conjecturedata]

Note: If you use TableForm as part of modeldata (at the end //TableForm), then the ListPlot will not work!

The graph looks linear. Of course, you know there is a “best” fit line for the data. Here we just guess the slope by  choosing two appropriate data points. The line through the far left data point (3929, 1379) and the one before the far right data point, seems to be a good approximation.  The slope is (32706.-1379)/(248710-3929) or approximately 0.128.

(32706.-1379)/(248710-3929)

Below is the way we will find the best fit line in future. Notice its slope is NOT very close to our estimated value.

Fit[conjecturedata,{1,x},x]

Our model is p[n+1]-p[n] = 0.128 p[n], or p[n+1] = 1.128 p[n]. Let’s see if it matches our data.

First define our function. Here Mathematica does not like defining a function with argument n+1. So we will use n and n-1.

Producing model data in this fashion will become impossible for more sophisticated models, as you will see in the next part!

p[0] = 3929;
p[n_] := 1.128  p[n - 1]
modelpop = Table[p[n], {n, 0, 22}]

modeldata=Transpose[{time,modelpop}]

modelplot=ListPlot[modeldata, PlotStyle->{Red}]

Show[modelplot, dataplot, PlotRange->All]

Our model plot, in red, does NOT seem to be a reasonable approximation.

If you plot both simultaneously, Mathematica will choose different colors for each one automatically. The first one in the list has the default (in this case navy) color.

ListPlot[{data,modeldata}, PlotRange->All]

For this homework problem, you should look for a better model! The above model implies that the line approximating the relationship between p[n+1]-p[n] and p[n] goes through the origin as indicated by p[n+1]=1.28 p[n]. But the Fit command suggested another line:   (p[n+1]-p[n]-4509) / p[n] = 0.099 or p[n+1] = 4509 + 1.099 p[n]. Is there a real-world justification for this?

Section 1.2, Example 2

time={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18}
pop={9.6,18.3,29,47.2,71.1,119.1,174.6,257.3,350.7,441,513,559,594.8,629.4,640.8,651.1,655.9,659.6,661.8}

Here I combine the data into ordered pairs.

data=Transpose[{time, pop}]

Just double checking that I have paired them correctly.

TableForm[data]

Below is the plot of the data.

dataplot=ListPlot[data]

pop[[n]] gives the nth data point in pop. Below I have formed the first differences.

popchange=Table[pop[[n+1]]-pop[[n]],{n,1,18}]

I will add a blank spot to the end of this list so it will be the same length as the other data and I use the transpose command to display all in table format.

popchange=Append[popchange,{}]

table=Transpose[{time,pop,popchange}]//TableForm

It looks like p[n+1]-p[n] is getting close to zero. Your book chooses the carrying capacity, largest population size, of 665. If growth was unbounded, we could say p[n+1]-p[n] is proportional to the p[n]; p[n+1]-p[n] = k p [n]. It is not the case here. Conjecture: Is p[n+1]-p[n] proportional to (665-p[n]) p[n];  p[n+1]-p[n] = k (665-p[n]) p[n]. If that is the case, the ratio (p[n+1]-p[n]) / (665-p[n]) p[n] should be the constant k. Plot (p[n+1]-p[n]) versus (665-p[n])p[n].

conjecturedata=Table[{(665-pop[[n]])pop[[n]],pop[[n+1]]-pop[[n]]},{n,1,18}]

Note: Pay attention to the order of data in the ordered pairs

TableForm[conjecturedata]

ListPlot[conjecturedata]

Note: If you use TableForm as part of modeldata (at the end //TableForm), then the list plot will not work!

The graph looks linear. Of course, you know there is a “best” fit line for the data. Here we just guess the slope by  choosing two appropriate data points. (0,0) should be a data point since if the original population is zero, then both terms should be zero. The far right data is the point (110225, 90.3). The slope is (90.3 - 0) / (110225 - 0)  or approximately 0.00082.

(90.3-0)/(110225-0)

Below is the way we will find the best fit line in future. Notice its slope is very close to our estimated value.

Fit[conjecturedata,{1,x},x]

Our model is p[n+1]-p[n] = 0.00082 (665-p[n]) p[n], or p[n+1] = p[n] + 0.00082 (665-p[n]) p[n]. Let’s see if it matches our data.

First define our function. Here Mathematica does not like defining a function with argument n+1. So we will use n and n-1. Although below is the normal way to define a function. For recursive equations, the time it takes to make the table will quickly become extremely long since for each n value, Mathematica will recalculate all the earlier ones! So, in the table below p[1] gets calculated 18 times, p[2] gets calculated 17 times, etc.

p[0] = 9.6;
p[n_] := p[n - 1] + 0.00082 (665 - p[n - 1]) p[n - 1]
modelpop = Table[p[n], {n, 0, 18}]

Instead, use RecurrenceTable since it remembers earlier values and won’t recalculate them!

(* Since p was used before, I "clear" it. You can use (* comments *) to add comments to your input cells, similar to adding comments to a program. *)
Clear[p]

modelpop=RecurrenceTable[{p[n+1]==p[n]+0.00082 (665-p[n]) p[n],p[0]==9.6},
p,{n,0,18}]

modeldata=Transpose[{time,modelpop}]

modelplot=ListPlot[modeldata, PlotStyle->{Red}]

Show[modelplot, dataplot, PlotRange->All]

Our model plot, in red, seems to be a reasonable approximation.

If you plot both simultaneously, Mathematica will choose different colors for each one automatically. The first one in the list has the default (in this case navy) color.

ListPlot[{data,modeldata}, PlotRange->All]

There is a way around time issue in using Table command in recurring relations, suggested by one of your classmates. His solution is to add the name to each level so the system keeps it in its memory and reuses it. In the following it is done by adding p[n] in the definition of the function p[n_].

p[n_] := p[n]= p[n - 1] + 0.00082 (665 - p[n - 1]) p[n - 1]

Run the following two cells to compare the time, in seconds, it takes for each one. The time is the first number in the output.  I have used maximum n value of 11, so you don’t have wait too long for the long version.

Clear[p];
p[0] = 9.6;
p[n_] := p[n - 1] + 0.00082 (665 - p[n - 1]) p[n - 1]
Timing[modelpop = Table[p[n], {n, 0, 11}]]

Clear[p];
p[0] = 9.6;
p[n_] := p[n]=p[n - 1] + 0.00082 (665 - p[n - 1]) p[n - 1]
Timing[modelpop = Table[p[n], {n, 0, 11}]]